AD-A190  568 


AFIT/GE/ENG/88M-3 


I 


MULTIPLE-MODEL  PARAMETER-ADAPTIVE  CONTROL 
FOR  IN-FLIGHT  SIMULATION 


THESIS 

Thomas  J.  Berens 
lLt.  USAF 


AFIT/GE/ENG/88M-3 


DTIC 


ELECTE 

MAR  2  8  1988 


/pprovoc  for  putlic  release ;  Distribution  ur limiter 


r> 

.  V 


A 


"  *V-  *V.  »  .  V.  ■ 


AFTT/GE/ENG/8  8M-3 


MULTIPLE  MODEL  PARAMETER  ADAPTIVE  CONTROL 
FOR  IN-FLIGHT  SIMULATION 


THESIS 


Presented  to  the  Faculty  of  the  School  of  Engineering 
of  the  Air  Force  Institute  of  Technology 
Air  University 

In  Partial  Fulfillment  of  the  Requirements  for  the  Degree  of 
Master  of  Science  in  Electrical  Engineering 


Thomas  j.  Berens,  B.S. 
First  Lieutenant,  USAF 


Accession  For _ 

ntis  gra&i  rjT 

DTIC  TAB  L> 

Unannounced  Q 

Justification- - ■ 

Bt - - - 

Distribution/  _ __ 

Availability  Codes 
lAvali  and/or 
Dlst  '  Special 


March  1988 


Approved  for  public  release;  distribution  unlimited 


V^*JT'.er<v'r'C'".C:  .  ^.rV^C^O'-CV'VrV'V  'Ti''’*  vV,rVaV <r,-y_  v.  V, v-.v,  rf\' y.  w-.'-'V  w  v.  y. y.  y_ «V  V.  vvww.r.K\ 


Preface 

The  purpose  of  this  study  was  to  modify  an  existing  parameter 
estimation  algorithm  to  a  multiple-input  multiple-output  difference 
equation  model.  The  immediate  need  for  this  procedure  was  for  use  in 
adaptive  control  of  an  in-flight  aircraft  simulator,  but  the  approach 
should  be  valid  for  other  applications. 

The  algorithm  was  used  in  conjunction  with  an  adaptive  control  law 
in  a  computer  simulation  of  an  in-flight  simulator  aircraft. 

Both  the  parameter  estimation  algorithm  and  the  control  law  worked 
well.  When  sensor  noise  was  added  to  the  simulation,  however,  estimator 
performance  suffered.  Future  study  must  modify  the  estimation  algorithm 
to  decrease  noise  sensitivity.  The  work  should  be  continued,  as  it 
provides  an  alternative  approach  to  fixed  gain  flight  control  systems 
which  require  gain  scheduling. 

This  work  was  a  direct  extension  of  an  earlier  thesis  by  Capt.  Luis 
A.  Pineiro  of  the  Flight  Dynamics  Tab  who  must  be  credited  with  the 
initial  concept  of  the  work.  I  am  deeply  indebted  to  Capt.  Pineiro  for 
many  hours  he  spent  explaining  various  aspects  of  parameter  identifica¬ 
tion  and  working  out  software  problems.  Finally,  I  would  like  to  thank 
my  advisor,  Col.  Daniel  Biezad  for  his  efforts  and  patience  during  a 
time  when  it  appeared  this  work  would  never  be  completed. 


Thomas  J.  Berens 
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Abstract 


Adaptive  control  of  aircraft  model-follcwing  systems  has  shown 
premising  results  for  in-flight  simulation,  but  the  computational 
expense  and  slew  convergence  of  conventional  parameter  estimation 
techniques  for  higher  order  models  inhibits  their  direct  use  for  in¬ 
flight  simulation.  Computer  simulations  of  adaptive  systems  usually 
assume  seme  knowledge  of  model  parameters  in  order  to  maintain  tracking 
fidelity  at  a  reasonable  computational  cost  as  parameters  change.  This 
thesis  incorporates  a-priori  information  into  a  multiple-model  estima¬ 
tion  algorithm  which  assigns  a  probability  weighting  of  each  estimator 
within  a  "bank"  of  estimators.  Final  parameter  estimates  used  in 
adaptive  control  are  formed  as  a  probabal istic  weighted  sum  of  individ¬ 
ual  estimates.  Simulations  of  the  system  shew  excellent  tracking 
performance  throughout  the  flight  envelope.  A  moving  bank  scheme  for 
use  over  a  wide  range  of  flight  conditions  is  recommended  as  a  further 
area  of  study. 


MJIUTPLE-MXEL  PARAMETER-ADAPTIVE  OOtOTOL  FOR  IN-FLIGHT  SIMJIATICN 

I.  Introduction 

An  in-flight  simulator  is  an  aircraft  whose  stability,  feel,  and 
flying  characteristics  can  be  changed  to  mimic  another  aircraft  from 
the  pilot's  perspective.  In-flight  simulators  are  used  in  new  aircraft 
system  pre-production  testing,  in  research  and  development,  and  in 
training  test  pilots  [3,  14]. 


Figure  (I— 1)  Vista  Concept 
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Figure  (1-2)  Flight  Envelope 


The  U.S.  Air  Force  currently  has  several  types  of  in-flight  simula¬ 
tors  in  its  inventory  and  will  soon  replace  its  current  fighter  simula¬ 
tor,  the  NT-33.  The  NT-33  was  developed  in  the  1950's  and  can  no 
longer  match  the  performance  of  new  fighter  aircraft.  Its  replacement, 
VISTA  (Variable  Stability  In-flight  Simulator  Test  Aircraft) ,  is  a 


— modern,  high  performance  fighter  aircraft,  modified  with  vari¬ 
able  stability  controls  and  a  reprogrammable  cockpit,  offering  the 
capability  to  perform  large  amplitude  maneuvers  over  a  wide  range 
of  operating  conditions  in  an  expanded  flight  envelope  [4] 


VISTA  has  spurred  renewed  interest  in  high  performance  adaptive 
control  since  it  is  an  ideal  testbed  for  new  flight  control  methods. 
Adaptive  control  is  well  suited  to  in-flight  simulation  due  to  the 


control  since  it  is  an  ideal  testbed  for  new  flight  control  methods. 
Adaptive  control  is  well  suited  to  in-flight  simulation  due  to  the 
time-varying  nature  of  the  host  aircraft's  stability  derivatives  as  the 
aircraft  transitions  within  its  flight  envelope  (see  Figure  1-2) .  This 
type  of  control  implies  that  estimates  of  the  stability  derivatives  are 
used  in  control  law  calculations  [21-23]  to  reduce  or  minimize  degraded 
tracking  performance  of  the  model  aircraft  by  the  host. 

Recursive  estimation  techniques  which  calculate  discrete  parameter 
values  from  input  and  output  samples  are  impractical  with  in-flight 
■  simulation  when  the  model  order  results  in  a  large  number  of  parameter 

i 

|  estimates  [25]  or  when  the  input  excitation  is  insufficient  [6]. 

I 

Pineiro  [19,20]  used  an  algorithm  developed  by  Hagglund  [7]  to  overcame 

the  excitation  problem,  but  limited  the  number  of  estimated  parameters 
i 

i  to  avoid  unacceptable  computational  expense  and  convergence  times.  In 

j  this  thesis,  known  information  about  the  parameters  is  combined  via  a 

multiple  model  estimation  algorithm  [1,2,15,26]  which  may  be  used  in 
conjunction  with  Hagglund' s  algorithm  to  overcome  these  limitations 
without  reducing  the  number  of  estimated  parameters.  Such  a-priori 
data  is  usually  available  from  flight  and  wind  tunnel  testing.  Use  of 

this  data  in  a  parameter-adaptive  model-follcwing  system  speeds  conver- 

i 

gence  and  significantly  reduces  the  computational  expense. 

Section  2  of  this  thesis  summarizes  the  model  following  system 
developed  by  Pineiro  which  implements  an  adaptive  fast  sampling  MEMO 
control  law  for  robust  tracking  and  which  partially  estimates  a  set  of 
linear  difference  equation  model  parameters.  Section  3  shews  hew  to 
estimate  all  parameters  at  a  relatively  small  computational  cost  by 
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using  a-priori  information.  The  multiple  model  estimation  algorithm  is 
emphasized  and  a  derivation  is  included.  Section  4  discusses  several 
modifications  to  the  multiple  model  algorithm  which  improve  stability, 
convergence  speed,  or  computational  cost.  In  Section  5  a  simulation 
demonstrates  the  estimation  of  parameters  which  are  blended  via  the 
multiple  model  algorithm,  filtered  [10],  then  fed  to  the  control  law. 
Finally,  conclusions  and  recommendations  for  further  study  follow  in 


Section  6. 
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I.  Parameter-Adaptive  Model-Following 

Figure  (II-l)  shews  a  block  diagram  of  a  parameter-adaptive  model- 
follcwing  system  [20] .  The  plant  is  a  MIMO  state  equation  which  simu¬ 
lates  the  longitudinal  dynamics  of  an  aircraft.  The  plant  is  control¬ 
lable  by  a  proportional -plus-integral  (PI)  control  law. 

This  section  describes  two  methods  of  calculating  the  gains  of  the 
control  law  [21-24].  The  fixed  gain  method  requires  knowledge  of  plant 
parameters  by  the  control  law.  When  these  parameters  are  unavailable, 
an  adaptive  method  is  used  to  estimate  them  from  real  time  input-output 
measurements  via  a  recursive  parameter  estimation  algorithm.  Several 
recursive  parameter  estimation  algorithms  are  presented  and  design 
issues  discussed. 


Figure  (II-l)  Adaptive  Model  Following  S inula ti on  Block  Diagram 


II. A.  Fixed  Gain  Control  Law 

The  aircraft  plant  is  represented  by  a  completely  controllable  and 
observable  MHO  state  space  model 

x(t)  =  Ax(t)  +  Bu(t) 
y(t)  =  cx(t) 
where 

x  (t)  e  Rml  a  e  r1^  y  (t)  e  r™*1 

u  (t)  e  R™*1  B  e  R™11  with  rank  "m" 

The  A,  B,  and  C  matrices  are  partitioned  according  to  the  control 
input  matrix,  B,  to  yield 


(2-la) 

(2-lb) 


where 

x^(t)  6  RPxl  B2  e  R(n-P) with  rank  "m" 

x2(t)  e  R(n_P)xl 


A  regular  plant  is  defined  as  having  a  first  Markov  parameter,  CB,  of 
full  rank.  Regular  plants  with  stable  transmission  zeroes  will  track 
input  given  the  discrete  output  feedback  control  law  (Figure  (II-2)). 


If  the*  matrices  A,  B,  and  C  in  eq.  (2-1)  are  known,  the  control 
gain  matrices  K]_  and  K2  in  eq.  (2-3)  are  calculated  by  discretizing  the 
plant  state  and  output  equations  in  block  diagonal  form  [20]  for  the 
sampling  period  T  as 


and 
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(2-6a) 


(2— 6b) 


(2— 6c) 


22  (k)  =  x2(k) 


(2— 6d) 


Ci  =  [  K!_1K2  ,  0  ] 


(2-6e) 


9 2  =  c2 


(2-6f) 
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ta12c2-1k1~1k2  >  In-m  +  A11T  "  TA12c2_lcl 


(2-6g) 
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The  input-output  relationship  can  be  expressed  in  terms  of  the 
closed-loop  transfer  function,  G(z) ,  where  y(z)  =  G(z)r(z)  and  z  is  the 
discrete  transform  operator .  As  the  sampling  frequency  is  increased 
G(Z)  assumes  the  asymptotic  form  [20] 


G (z)  =  G!(z)  +  G2(z) 

where  z  is  the  discrete  transform  operator  and 


Gf(z)  - 
G2(z)  = 


with 


C!(zln  -  In  -  TAq) ~ 1TB0 
^2(zIm  “  -^m  ~  ^4)  ^^4 


Kr1K2 

Al2C2_1Kr1K2  , 


a11_a12c2-1c1 


(2-7a) 


(2— 7b) 
(2-7c) 


(2-7d) 


Bq  = 


a12C2' 


A4  =  -b2k1c2 

b4  =  B2Kx 


{2-1 e) 


(2-7f) 

(2-7g) 


G^Cz)  and  G2(z)  are  the  slow  and  fast  mode  transfer  functions 
respectively.  The  slow  modes  can  be  grouped  into  two  sets  Z-y  and  Z2 
and  are  given  by 


Zx  =  {z€C:  det{ zlm  -  Im  +  TKi'1^}  =  0}  (2-8) 

z2  =  {z  €  C:  det{zln_m  -  In_m  -  TAn  +  TA12c2_lcl>  =  °>  (2-9) 

The  fast  modes  are  given  by 


Z3  =  {z  6  C:  det{zlm  -  Im  +  C2B2K1}  =  0} 


(2-10) 


Because  of  the  form  of  AO,  BO,  and  Cl,  the  eigenvalues  of  AO  are 
uncontrollable  or  unobservable.  Thus,  as  the  sampling  frequency 
increases,  the  slow  transfer  function  asymptotically  approaches  zero 
and  the  overall  system  transfer  function  contains  only  the  fast  modes, 
as  given  by  G2 (z)  which  can  be  put  in  the  equivalent  form 


G(z)  -  G2 (z)  -  (zlm  -  Im  +  C2B2K1) -1C2B2Kl 


The  controller  matrices  K]_  and  K2  are  then  given  by 


(2-11) 


Kl  =  [C2B2]“1S 
K2  =  qKx 


(2-12a) 

(2-12b) 


where  q  is  any  positive  scalar  greater  than  zero,  and  S  is  a  diagonal 
tuning  matrix.  Both  q  and  S  are  chosen  by  the  designer  to  achieve  the 
desired  tracking  characteristics. 


II. B.  Adaptive  Control  Law 

If  A,  B,  and  C  in  eq.  (2-1)  are  unknown,  the  control  law  gain 
matrices  of  eqs.  (2-3)  can  be  calculated  adaptively  by  expressing  eq. 
(2-8)  as 


y(k)  =  Ad(T)y(k  -  1)  +  H(T) u(k  -  1) 


(2-13) 


where  Ai(T)  =  exp{AT)  and  H(T)  is  defined  as  the  step  response  matrix. 
H(T)  can  be  calculated  in  terms  of  the  continuous  time  state  space 
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system  as: 


fT 

H(T)  =  I  C  exp(At)  B  dt  (2-14) 

J  0 

* 

Expanding  exp (At)  as  a  Taylor  series  about  t  =  0  and  dropping  the  high 
order  terms  yields  for  small  sampling  periods  H(T)  =  TCB,  and  the 
control  gains  can  be  expressed  as 

Kx  =  H"1(T)S  (2-17a) 

K2  =  qKx  (2— 17b) 

The  step  response  matrix  H(T)  can  be  estimated  from  real-time  input- 
output  data  by  expressing  eq.  (2-13)  in  terms  of  an  nth  order 
autoregressive  vector  difference  equation  of  the  form 


y(k)  =  Bxu(k  -  1)  -  Axy(k  -  1)  +...+  BnU(k  -  n) 

-  Afjy(k  -  n)  +  e(k)  (2-16) 

where  the  equation  error  vector  e(k)  is  assumed  to  be  a  zero-mean 
Gaussian  white-noise  vector  with  variance  v(k)  and  the  matrices  A^  e 
fjmxm  (i=if  2,  —  ,n),  B^  e  r1™11  (i=l,  2,...,n)  are  the  parameters  of 
the  nth  order  difference  equation. 

Eq.  (2-16)  can  be  expressed  alternatively  as 


y(k)  =  ©oT(k)  +  e(k) 


(2-17a) 


where,  for  a  SISO  system, 


0T(k)  =  [-yT(k  -  l),...,-yT(k  -  n),uT(k  -  l),...,uT(k  -  n) ] 

e  R1*2  (n)  (m)  (2-17b) 

and 

eT(k)  =  ...,Bn]  eR^WW  (2-i7c) 

The  step  response  matrix  is  updated  by  invoking  the  certainty 
equivalence  principle  [20].  By  definition  of  the  step-response  matrix 
it  can  be  shown  that  [22,23] 

H(T)  =  TCB  =  B1  (2-18) 

All  the  parameters  of  eq.  (2-16)  must  be  estimated  to  identify  Bj_. 

The  input-output  relationship  of  eq.  (2-6)  is  transformed  into  eq.  (2- 
16)  via  the  following  steps: 

Step  1:  Obtain  the  Z-transform  of  (2-6) 

zx(z)  =  Acjx(z)  +  E^ju(z)  (2-19a] 

y(z)  =  Cx(z)  (2-19b) 

Step  2:  Obtain  an  input-output  relationship  independent  of  x(z): 


(zl  -  Ad)x(z)  =  BdU(z) 

(2-20) 

x(z)  =  (zl  -  Ad)-1E^(z)u(z) 

(2-21) 

y(z)  =  G(z)u(z) 

(2-22) 

12 


where 


G(z)  =  C(zl  -  Ad)-lEti(z) 


G(z)  is  a  polynomial  matrix  of  the  form: 


Gn (z)  G^2(z)  •••  G]j^(z) 

G(z)  =  G2i(z)  G22(z)  —  G2m(z) 

•  •  t 

•  •  * 

Gml(z)  Gj^2(z)  ...  Gjfl^z) 


(2-23) 


(2-24) 


where  G^j (z)  is  the  transfer  function  relating  the  output  y^(z)  to  the 
control  input  uj(z)  and  is  of  the  form 


fc>izw  +  b2Zw  1  +  ...  +  byZ  + 


Gij(z)  = 


zn  +  axzn  1  + 


+  an-lz  +  an 


(w  <  n)  (2-25) 


by  dividing  each  numerator  and  denominator  in  G(z)  by  zn,  G(z)  is 
transformed  into  a  delay  operator  of  the  form 

b^zw-n  +  b2Zw-n_1  +  —  +  b^+iz-11 

Gij(z)= -  (w  <  n)  (2-26) 

1  +  a^z-1  +  . . .  +  an_iz  +  anz  n 


Step  3:  Redefine  the  input-output  relationship  obtained  in  step  2, 
expressed  as  eq.  (2-22),  in  terms  of  polynomial  matricies,  A(z)  and 
B(Z) ,  i.e.: 


A(z)y(z)=B(z)u(z) 


where 


A(z)=  M  +  A^z-1  +  ...  +  Anz  n 


(2-27a) 


(2-27b) 
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B(z)  =  M  +  +  ...  +  B^+iZ-1"1 


(2-27 c 


Aii  - 


aii{ll) (z)  aii( 12 } (z)  ...  aii(lm}(z) 
aii{21> (2)  aii(22}(z)  ...  aii(2m>  (z) 

aii(ml)(z)  aii(m2}(z)  ...  ^{mm}  (z) 


(2-27d; 


where  is  the  a-[i  coefficient  of  Gij(z),  eq.  (2-26)  and  M  is  an 

m  x  m  matrix  of  one's.  Bj_j_  are  of  the  same  form  as  eq.  (2-24d) . 

Step  4:  Obtain  the  inverse  Z-transform  of  eq.  (2-27)  and  rearrange 
terms  to  obtain  the  linear  difference  equation  model  given  by  eq.  (2- 
16)  . 


II. C.  Recursive  Parameter  Estimation 

The  parameter  vector,  9,  of  eq.  (2-17)  can  be  obtained  from 
input-output  measurements  using  a  recursive  parameter  estimation 
algorithm.  The  recursive  least  squares  (RLS)  algorithm  is  a  basic, 
widely  used  parameter  estimation  method  [7,13].  The  RLS  algorithm 
assumes  0  is  a  random  vector  with  a  Gaussian  prior  distribution  of 
mean=9(k  -  1)  and  variance  P(k  -  1) .  Input-output  measurements  are 
correlated  with  ©,  so  at  time  kT  the  posterior  probability  density 
function,  p(©|yk,  uk) ,  can  be  formed,  where 

yk  =  (y(k) ,  y(k  -  l) , . . . ,y(i) ) 
uk  =  (u(k) ,  u(k  -  1) , . . . ,u(l) ) 

9  can  be  recursively  updated  for  a  SISO  system  as 

9(k)  =  9(k  -  1)  +  (l/v(k))P(k)o(k)e(k)  (2-28) 


P(k  -  l)0(k)oT(k)P(k  -  l) 

P(k)  =  P(k  -  1)  -  -  (2-29) 

v(k)  +  oT(k)P(k  -  l)0(k) 

where 

v(k)  =  variance  of  e(k) 

P(k)  =  estimated  parameter  covariance  matrix  e 
The  proof  [13]  is  based  on  Bayes  rule  in  the  form 

P(A  j  B,  C)  =  P(B|A,C)P(A|C)/P(B|C)  (2-30) 

Where  P(A|B,C)  is  the  probability  of  the  event  A,  conditioned  on  B  and 
C.  Assuming  is  deterministic  and  applying  this  formula  to  the 
posterior  density  gives 

P(©lYk)  =  p(0|y(k) ,yk_1) 

=  P(y(k) |0,yk“1)p(©|yk-1)/p(y(k) |yk-i)  (2-3i) 

The  desired  result  is  proved  by  induction. 

Step  1; 

p(0|yo)=(27T)-dim  6/2(det  P(0) ) ~1/2exp{ _i/2 [e 

-  6(0) ]Tp_1[0  -  0(0) ] }  (2-32) 


by  assunption. 

Step  2:  Assume  that 

P(©l  yk_1)  =  (2  7T)'<3iim  e/2(det  P(k  -  1)  )_1/2exp{-l/2[0 

-  6(k  -  1) ]TP-1 (k  -  1) [6  -  0(k  -  1) ] }  (2-33) 


Now  calculate  p^ly1-)  using  eq.  (2-31).  From  eg.  (2-17) 


Under  the  assumption  that  e^  is  a  sequence  of  independent  random 
variables  with  zero  means  and  variances  v(k) , 

p(y(k)  |e,yk-1)  =  (27Tv(k)  )-1/2e>p{-[l/2v(k)  ]  [y(k) 

-eTo(k)]2}  (2-35) 

Hence,  from  eq.  (2-31) , 

p(9|yk)  =  Norm  exp{-[l/2v(k) ] [y(k)  -eTo(k)]2 
-1/2 [©  -  0(k  -  1) ]TP_1(k  -  l) [0  -  e(k  -  1) ] )  (2-36) 

where  the  0  independent  normalization  factor  is  not  explicitly 
written  out  . 

The  exponent  is  now  written  as  a  quadratic  form: 

-2  log  p(0|yk)  =  const  +  [l/v(k) ]y2 (k) 

-  [l/v(k)]y(k)©T(k)0  -  [l/v(k)]©T0(k)y(k) 

+  [ l/v (k) ] 0To ( k) oT (k) ©  +  eTP-1 (k  -  1)0 

-  ©TP-1(k  -  l)9(k  -  1)  -  0T(k  -  l)P_1(k 

-  1)9  +  ©T(k  -  l)P_1(k  -  l)9(k  -  1)  (2-37) 

Define 

P_1  =  P_1(k  -  1)  +  [l/v(k) ]o(k)0T(k) ;  (2-38) 

thus  the  preceding  expression  is  written  as 

const  +  [l/v(k) ]y2 (k)  +  ©(k  -  l)Tp-l(k  -  l)©(k  -  1) 

+  0TP_1©  -  ©T[ [ 1/v (k) ]o(k)y(k)  +  P-1(k  -  l)9(k  -  1)] 


-  [[l/v(k)]e>(k)y(k)+p-l(k  -  l)0(k  -  1)]T© 

=  const'  +  [0  -  P[l/v(k)  ]o(k)y(k)  -  PP-1(k  -  l)©(k  -  1)  ]T 
P-1]]©  -  P[l/v(k) ]0(k)y(k)  -  PP-1(k  -  l)©(k  -  1)] 


where  const '  is  a  new,  ©  independent  normalization  constant.  Since 


PP-1 (k  -  1)  =  I  -  [l/v(k) ]P0(k)oT(k) , 


The  expression  within  the  parentheses  nay  be  written  as  ©  -Q, 
where 


e  =  ©(k  -  1)  -  [ l/v(k,  ;PoT(k) ;y(k)  -  ©T(k  -  l)o(k)]. 


This  means  that  the  posterior  density  at  time  kT,  p(©|yt),  is 
Gaussian  with  mean^©,  given  by  eq.  (2-41),  and  covariance  matrix  P, 
given  by  eq.  (2-42) .  Applying  the  matrix  inversion  lemma*  to 
eq.  (2-35)  yields 


P  =  P(k  -  1)  -  P(k  -  l)0(k)oT(k)P(k  -  l)/[v(k) 
+  oT(k) P(k  -  l)o(k) ] 


(2-42 


Eq.  (2-28)  and  eq.  (2-29)  can  also  be  derived  by  minimization  of  the 


The  Matrix  Inversion  lemma  states  that 
[A  +  BCD]-1  =  A-1  -  A-1B[DA-1B  +  C-1]-1DA" 


where  A,  B,  C,  and  D  are  matrices  of  compatible  dimensions,  so  that  the  product 
BCD  and  the  sum  A  +  BCD  exists. 


cost  function  [13,20]. 


J{9(k)}  =  SUM{v-2(j)[y(j)  -0(j)0(j)]2}  (2-43) 

j=l 


II. D.  Estimating  v(k] 

The  variance  of  the  prediction  error  is  not  explicitly  updated  by 
eq.  (2-28)  or  eq.  (2-29) .  Given  the  probability  distribution  of  the 
prediction  error  is  white,  then  an  estimate  of  the  variance  can  be  made 
from  the  sequence  e^-1  based  on  the  central  limit  theorem: 


v(k)  =  [l/k]SUM{e(j)2} 

j=l 


(2-44) 


If  v(k)  is  slowly  time  varying,  then  v(k)  can  be  tracked  with  a  fading 
memory  filter  of  the  form 


v(k)  =  mu{v(k  -!)}+(!-  mu)e2(k) 


where  0  <  mu  <  1 


(2-45) 


Ihe  smaller  the  value  of  mu,  the  more  information  is  forgotten  during 
each  update. 


II. E.  Modifications  to  KLS  for  In-flight  Simulation 


Ihe  RLS  algorithm  given  by  eq.  (2-28)  and  (2-29)  assumes  constant 


parameters  which  are  actually  a  function  of  an  aircraft's  flight 


condition.  A  simple  method  for  tracking  time-varying  parameters  is  to 
exponentially  discard  old  information  by  incorporating  a  constant 
forgetting  factor  into  eq.  (2-29)  [7,13],  i.e. 


P(k  -  l)0(k)oT(k)P(k  -  1) 

P(k)  =  (l/g)P(k  -  1)  -  -  (2-46) 

(l/g)v(k)  +  oT(k)P(k  -  l)0(k) 


where  0  <  g  <  1.  Choosing  a  value  of  g  is  a  trade-off  between  tracking 
time  varying  parameters  and  discarding  good  information.  A  small  value 
of  g  allows  the  R1S  algorithm  to  quickly  track  changing  parameters. 

Accurate  parameter  estimation,  however ,  requires  sufficient  input 
excitation  [6],  If  current  measurements  contain  no  new  useful 
information  and  old  information  is  being  discarded,  parameter  estimates 
may  "burst"  from  the  best  fitting  parameter  values.  This  phenomena  is 
referred  to  in  the  literature  as  estimator  wind-up  [27]. 

Hagglurd's  algorithm  [7]  attempts  to  solve  the  wind-up  problem  by 
-.mounting  past  information  in  such  a  way  that  if  the  parameters  were 
constant,  a  constant  amount  of  information  is  retained  [7]. 
Specifically,  Hagglurd's  algorithm  modifies  the  P(k)  update  equation 
as 


P(k)  =  F(k  -  1) 

P(k  -  l)o(k)0T(k)P(k  -  1) 

[v“i (k)  -  a(k) ]-1  +  0T(k) P(k  -  l)o(k) 


(2-47) 


where  a(k)  is  a  discounting  factor  dependent  upon  the  amount  of 
excitation  present  in  the  input  and  output  signal.  The  entire 
algorithm  is  listed  in  Appendix  C. 


The  RLS  algorithm  with  Hagglund's  modifications  to  eq.  (2-28) 
converges  slowly  for  high  order  models  requiring  a  large  number  of 
parameter  estimates.  Fault  detection  is  a  method  which  may  decrease 
convergence  time.  A  fault  is  defined  as  a  rapid  change  in  system 
dynamics.  After  a  fault  the  current  parameter  estimates  no  longer 
accurately  describe  aircraft  dynamics.  The  estimation  algorithm 
converges  but  only  after  an  undesirable  time  delay  which  degrades 
control  performance. 

Hagglund's  algorithm  accounts  for  the  loss  of  parameter  knowledge 
by  increasing  the  value  of  the  estimated  parameter  error  covariance 
when  a  fault  is  detected.  The  increased  weight  placed  on  new 
measurements  allows  the  estimation  algorithm  to  converge  faster.  A 
fault  is  expressed  by  modifying  eq.  (2-46)  as 


P(k)  =  P(k  -  1) 

P(k  -  l)o(k)0T(k)P(k  -  1) 

-  +  B(k)  (2-48 

[v"1(k)-a(k) ]-1  +  0T(k)P(k  -  l)o(k) 

where  B(k)  is  a  matrix  which  depends  upon  the  size  of  the  fault  which 
increases  P(k)  after  a  large  parameter  change  occurs.  Criterion  for 
determining  when  a  fault  has  occurred  and  how  to  choose  B(k)  is 
discussed  by  Hagglund  [ 7 ] . 

Hagglund's  algorithm  alone  is  still  inadequate  for  many  higher 
order  systems.  Therefore,  Pineiro  assumed  elements  of  the  parameter 
vector  9(k)  not  used  in  the  control  law  remain  constant  about  a 
nominal  flight  condition.  In  such  case  these  elements  are  not 


estimated  and  eq.  (2-17)  must  be  partitioned  as 


(2-50a) 
(2— 50b) 
(2— 50c) 


y (k)  =  eTv(k)0v(k)  +  ©TcOcCk)  +  e(k) 

©T(k)  =  [cvT(k)  ;  0cT(k) ] 

eT(k)  =  [0vT(k)  ?  0cT(k) ] 

where  the  subscript  v  denotes  the  identified  parameters  and  c  denotes 
the  constant  parameters.  Note  0c(k)  corresponds  to  9C  but  is  not 
constant.  The  relationship  between  6v(k)  and  B]_  is  illustrated  in 
the  following  example. 

Example  2-1 

Assume  a  linear  difference  equation  of  the  form 

y(k)=  -Ajylk  -  1)  -  A2y(k-2)  +  Bxu(k  -  1)  +  B2u(k-2)  +  e(k)  (2-51) 

where  y(k)  e  R2xl,  Ax  6  R2x2,  A2  6  R2x2 

u(k)  e  r2x1,  Bx  e  r2x2,  B2  e  r2x2 

This  equation  can  be  partitioned  in  the  form  of  eq.  (2-51)  to  isolate 
Bi  for  identification,  i.e. 


®cT(k) 

=  [-yT 

(k  - 

1), 

-yT(k-2),  uT(k-2)]  6  Rlx6 

(2-52a) 

*cT(k) 

=  uT(k 

-  1) 

e  r1x2 

(2 -52b) 

ecT(k) 

=  [Ai, 

a2. 

B2] 

e  Rlx2 

(2-52 c) 

©vT(k) 

=  BX 

(2-52d) 
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This  section  discussed  several  basic  and  important  design  issues  in 
simulated  adaptive  control  of  an  in-flight  simulator.  The  section 
began  by  describing  the  plant  and  its  control  by  a  fixed  gain  PI  con¬ 
trol  law.  This  simulation  was  then  modified  by  replacing  the  fixed 
gain  controller  with  an  adaptive  system.  An  adaptive  system  bases 
control  gain  calculations  on  input-output  measurements  rather  than 
explicit  knowledge  of  system  parameters.  Adaptive  systems  are  more 
complex  since  they  require  a  parameter  estimation  algorithm.  A  well 
known  and  widely  used  recursive  parameter  estimation  algorithm,  RLS, 
was  derived.  Modifications  to  the  RLS  algorithm,  however,  are  required 
for  use  in  in-flight  simulation.  Section  3  extends  recursive  parameter 
estimation  to  make  better  use  of  available  a-priori  information. 


Parameter  Estimation  Us. 


A-Priori  Information 


im 

Standard  recursive  parameter  estimation  algorithms  are  not  well 
suited  to  adaptive  flight  control.  The  computational  expense  and  slow 
convergence  of  the  RLS  algorithm  makes  estimating  all  the  parameters  of 
higher-order  linear  difference  equation  aircraft  dynamics  models  which 
require  a  large  number  of  parameter  estimates  impractical  for  on-line 
use.  Even  modifying  standard  recursive  parameter  estimation  techniques 
by  adding  a  fault  detection  algorithm  are  found  to  be  inadequate  with¬ 
out  assuming  partial  knowledge  of  model  parameters  [20]. 

This  section  develops  the  multiple  model  algorithm  (MMA)  which  is 
well  suited  to  adaptive  flight  control.  It  estimates  parameters  from  a 
finite  amount  of  a-priori  information  which  speeds  parameter  estimator 
convergence  and  reduces  computational  expense.  A-priori  information 
available  from  wind  tunnel  and  flight  testing  is  in  the  form  of  models 
of  aircraft  at  various  flight  conditions.  The  parameter  vectors  of  eq. 
(2-16)  are  calculated  for  each  model  using  the  procedure  described  in 
Section  2  and  are  incorporated  into  the  MMA  as  means  of  Gaussian 
distributions.  The  final  model  parameter  vector  used  by  the  control 
law  is  a  random  variable  whose  conditional  distribution  function 
(conditioned  on  output  prediction  errors)  is  a  weighed  sum  of  these 
Gaussian  distributions.  Using  Bayes  rule,  the  weighting  coefficients 
are  updated  recursively  as  functions  of  output  prediction  errors  and 
estimates  of  the  output  prediction  error  variance. 


% 


III. A.  Finite  Set  Parameter  Estimation 

Before  discussing  the  MMA,  consider  the  more  general  case  of 
parameter  estimation  from  a  finite  set  of  models.  Assume  9(k)  is 
limited  to  nm  "candidate"  model  parameter  vectors,  ©i(k) ,  such  that 

nm 

e(k)  =  suMjaiCkjeiCk) >  (3-1) 

i=l 

where  ©i(k)  represents  the  parameter  vectors  of  linearized  input- 
output  models  of  a  non-linear  systems  at  nominal  operating  points. 

Gain  scheduling,  a  method  often  used  to  choose  a  dynamics  model  for 
flight  control,  is  an  example  of  parameter  estimation  from  a  set  of 
candidate  models.  Gain  scheduling  chooses  aircraft  dynamics  models 
based  on  dynamic  pressure  which  is  a  function  of  an  aircraft's  airspeed 
and  altitude.  A  disadvantage  of  gain  scheduleing  is  dynamic  pressure 
is  measured  with  external  sensors. 

Many  quantities  are  available  which  require  no  external  sensors 
such  as: 

1)  Ihe  parameter  covariance  matrix  (3-2) 

Pi(k)  =E<(e(k)  -9(k))(0(k)  -0(k))T) 

2)  The  variance  of  the  prediction  error  (3-3) 

vf(k)  =  E{ei(k)eiT(k) ) 

3)  Ihe  Autocorrelation  Function  (3-4) 

Aee(k,j)  =  E{ei(k)eiT(k  -  j) } 

4)  Ihe  Crosscorrelation  Function 
Aue(k,j)  =  E{u(k)eiT(k  -  j)  ) 


(3-5) 


III.B. 


Comparison  of  Parameter  Fit  Quantities 
If  each  parameter  vector  of  eq.  (3-1)  is  estimated  with  a  recursive 
parameter  estimation  routine  such  as  RLS,  the  parameter  variance  can  be 
used  as  a  model  fit  indicator.  Unfortunately,  basing  a  model  selec¬ 
tion  test  on  the  parameter  variance  limits  the  designer.  For  example, 
when  model  parameters  vary  slowly  with  time,  Pi(k)  may  be  modified  non- 
linearly,  as  discussed  in  Chapter  2,  to  prevent  estimator  wind-up  and 
optimize  parameter  tracking  performance.  Another  example  occurs  in 
fault  detection  schemes.  Parameter  estimates  converge  faster  after  a 
large  parameter  change  when  P^(k)  is  enlarged.  Both  these  schemes 
destroy  goodness  of  fit  information  contained  within  P|(k) ,  since  eq. 
(3-2)  no  longer  holds  true. 

Prediction  errors,  e^(k),  are  not  effected  seriously  by  ad-hoc 
parameter  tracking  schemes.  Also,  a  parameter  estimation  algorithm 
need  not  be  used  to  calculate  0(k)  .  This  greatly  reduces  the  compu¬ 
tational  cost  involved  in  calculating  9(k).  Estimates  of  the  predic¬ 
tion  error  variance,  however,  may  be  masked  by  noise. 

A.  The  Multiple  Model  Algorithm  (MMA)  -  Derivation 
Ihe  multiple  model  algorithm  (MMA)  provides  a  method  of  blending  a- 
priori  data  with  parameter  estimation.  Although  similar  to  the  RLS 
algorithm,  this  approach  views  9(k)  as  a  random  variable  whose 
probability  density  function  (pdf)  is  approximated  as  a  weighted  sum  of 
nm  Gaussian  pdf's  rather  than  a  single  Gaussian  pdf. 
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yyA'.o 


^ 


p(9)  =  (27T)“dim  e/2-V2  SUM{ai(0)  [det(Pi(0))  ]_1/2 

i=l 

exp{-l/2[0  -  0i(O)]TPi-l(O)[©  -  6i(0)]} 


where  ©i(0)  are  parameter  vectors  such  as  those  in  eq.  (3-1) ,  which 
represent  the  means  of  each  Gaussian  pdf,  and  Pi(k)  represents  the 
variance.  The  ai(0)  are  weighting  coefficients  where  0  <  ai(0)  <  1 


SUM{ai(0) )  =  1 
i=l 


Eq.  (3-6)  can  be  recursively  updated  from  input-output  data  via  the 
multiple  model  algorithm  where: 

m 

P(©lyk~1)  =  (2  7T)_dim  ©/2-1/2  sUM{ai(k  -  l)[det(Pj_(k  -  l))]*1/2 

i=l 


exp{-l/2[©  -  ©i(k  -  l)]TPi-l(k  -  1) [© 
-  ©i(k  -  1)]) 


p(e|yk)  =  (27T)_diln  ©/2-1/2  sUM{ai(k)/det(Pi(k) )  ]_1/2 

i=l 


exp{ -1/2 [©  -  ©i(k) ]TPi"1(k) [©  -  ©i (k) ] } 


where 


Pi(k)  =  Pi(k  -  1)  -  Pi (k  -  l)o(k)[v(k) 

+  o(k)TPi(k)ofk) ]_1o(k)TPi(k  -  1) 

©i(k)  =  ©i(k  -  1)  +  [1  v(k) ][Pi(k)o(k)ei(k) ] 
ai (k)  =  Csai(k  -  l)[v(k)  +  0(k)TPi(k)0(k)]-V2 

exp{-l/2eiT(k)  [v(k)  +  o(k)TPi  (k)o(k)  ]-:ie,  (k) 


where  Cg  is  a  normalization  factor,  such  that 


nm 

SUM{aj.(k) )  =  1 
i=l 

(3-13) 

nm 

0(k)  =  SUM(ai(k)0i(k) 

(3-14) 

i=l 

>i(k)  -  y(k)  -  ©iT(k)o(k) 

(3-15) 

The  proof  (Andersson)  which  is  similar  to  the  RLS  proof  in  Section  2  is 
based  on  Bayes  rule.  Applying  this  formula  to  the  posterior  density 
gives 

p(9|yk)  =  p(0|y(k) ,yk_1) 

=  p(y(k) |e,yk-1)p(0|yk-1)/p(y(k) lyk_1)  (3-i6) 

Now  calculate  p(0 | yk)  using  eq.  (3-16).  From  eq.  (2-19)  obtain 

e(k)  =  y(k)  -  0To(k)  (3-17) 

Under  the  assumption  that  (e(k) )  is  a  sequence  of  independent  random 
variables  with  zero  means  and  variances  v(k) ,  obtain 

P(y(k) |e,yk_1)  =  (27rv(k))_1/2exp{-[i/2v(k) ] [y (k) 

-0To(k)]2)  (3-18) 

Combining  eq.  (3-8)  and  eq.  (3-18)  via  Bayes  rule  yields: 


P(©lY*) 

nm 

=  C[v(k)  ]~V2  (2  7T)  (9/2  }-l/2  SUM{  ai  (k)  [det  (Pj_  (k  -  l))]-1/2 

i=l 

exp{ -1/2 [6  -  ©i(k  -  1) ]TPi-1(k) [0  -  9i(k  -  1)] 
exp{-[l/2v(k)][y(k)  -  eT0(k)]2)  (3-1S 

Equation  (3-9)  is  new  equated  with  eq.  (3-19)  to  find  matrices, 
0|(k)  and  Pi(k),  and  weighting  factors,  a^(k),  of  eq.  (3-9).  For 
notational  convienience ,  the  following  substitutions  are  made:  Pj_ 
Pi(k-l);  P'i  =  Pi(k);  e'i  =  6i(k);  ©i  =  6i(k  -  1);  a'i  = 
a^(k) ;  a^  =  a^(k  -  1) . 

The  exponent  terms  of  eq.  (3-9)  is: 

[e  -  e'jjTp'i-1!;©  -  e'jj  =  eTP'i-I9  - 

29' iTp' i_10  +  eTP'i_1e  (3-2 


and  the  exponent  terms  of  eq.  (3-19)  is: 

[ i/v ] [ y  -  eT0]2  +  [9  -  0i]TP’i-1[9  -  ©i] 

=  [y2  -  2y©Te>  +  (0T0)2][l/v]  +  0Tpi-1©  -  20iTpi'19 
+  0iTPi_19i  (3-21 

inspection  of  eq.  (3-20)  and  eq.  (3-21)  new  gives: 


[ -2y®/v  -  20jPj-1]0  =  -2e'iTP'i_19 


(3-2 


eT[0OT/v  +  Pi-1]©  =  ©Tp'i-1©  (3-23) 

eq.  (3-22)  and  eq.  (3-23)  must  be  satisfied  for  all  ©,  hence: 

©'i  =  P'i[Pi_1«i  *-  (l/v)oy]  (3-24) 

pi'_1  =  pi_1  +  ( 1/v)  e>0T  (3-25) 

Replacing  Pi-1  in  eq.  (3-24)  by  eq.  (3-25)  gives: 

9 ' i  =  ©i  +  (l/v)P'i-1oe  (3-26) 

where,  e  =  y  -  ©To 

Applying  the  matrix  inversion  lemma  to  eq.  (3-25)  gives: 

p,i  =  pi  -  PiOOTPi/[v  +  0TP-i0]  (3-27) 

Inserting  eq.  (3-23)  and  eq.  (3-24)  into  eq.  (3-20)  and  also  using  eq. 
(3-27)  gives 

[©  -  ©i]TP'i_1[©  -  ©i]  =  ©TPi-1©  -  2©iTPi-1©  + 

(1/V) (©To)2  -  (2/v) oT0y  +  ©iTPi_1©i  +  (l/v)(2©iT  + 

(l/v)oTPiy)oy  -  [l/(v  +  oTPi®) ] [ (©T0) 2 

+  2 ( 1/v) ©iTooTPioy  +  (1/v) (oTPj©y)2]  (3-28) 


The  first  five  terms  of  the  right  hand  side  of  eq.  (3-28)  are  identical 
to  eq.  (3-21),  except  for  the  term  y2(l/v)  in  eq.  (3-21).  Subtracting 


the  right  hand  side  of  eq.  (3-28)  for  the  right  hand  side  of  eq.  (3-21) 
new  gives  the  remaining  part: 


(l/v)y2  -  (1/v) (2©iT  +  (l/v)oTPiy]oy  +  [v  +  ©TPiy©]-1{ (©To) 2  + 
2 ( 1/v) e^TooTPi©y  +  ( l/v) 2 (oTPj©y) 2 ) 

=  (y  -  ©iTo)2/[v  +  0TP^o]  =  e^2/[v  +  ©TPi©] 

Eq.  (3-19)  can  new  be  written: 

m 

Cv~1/2(2  )-dim(©)/2-l/2  sUM{a;j_exp{ -(1/2) e-j2 [v  +  ©Tpj©) 

i=l 

det{ P^-1/2 )exp( -1/2 [©  -  ©'i]TP'i-1[©  -  ©'i]} 

Canparing  with  eq.  (3-9)  gives: 

a'i  =  Cai[det{P'i)/det{Pi)]V2e>p{-l/2ei2[v  +  ©TPi©] } 

Using 

detfP'i)  =  detjPi  -  Pj©0TPi/[v  +  ©Tpj©] } 

=  det{Pj_)[l  -  ©TPj0/[v  +  ©TPi©]  ] 

gives 

a'i  =  Cai[v/[v  +  ©TP^©]  ]  V2exp{ -(l/2)ej2/[v  4  oTPj_©]  } 

Hence,  a'i  cam  be  obtained  by  first  computing 


(3-29) 


(3-30) 


(3-31) 


(3-32) 


2i  =  ai[v  +  oTPiO]_1/2exp{-l/2ei2/[v  +  0TPi©] } 


(3-33) 


for  i  =  1, _ ,m  and  then  setting 

ran 

a'i  =  [SUM  aj]_1/2ai  (3-34) 

j=l 

This  section  developed  the  MMA  as  a  method  of  incorporating  avail¬ 
able  a-priori  information  into  an  adaptive  control  process.  The  MMA 
was  shown  to  be  similar  to  a  multiple  number  of  RLS  estimators  running 
in  parallel.  This  observation  is  used  in  the  following  section  to 
develop  several  mechanizations  of  the  MMA  used  for  simulated  on-line 
adaptive  flight  control. 


IV. A.  Second; 


arv  Estimator  Structure 
Comparison  of  eq.  (2-28)  and  eq.  (2-29)  with  eq.  (3-10)  and  eq.  (3- 
11)  shows  the  structure  of  the  secondary  parameter  estimate  (©f(k)) 
update  is  that  of  the  RLS  algorithm.  Therefore,  eqs.  (3-10)  and  (3-11) 
can  also  be  substituted  with  Hagglund's  Algorithm.  The  number  of 
parameters  which  need  be  updated  via  a  recursive  parameter  estimation 
routine  is  another  design  option.  A  comparison  of  estimating  all,  some 
(as  in  Pineiro's  simulation) ,  or  none  of  the  secondary  parameters  will 
be  made. 


TV. A. 1.  Full-Scale  Secondary  Parameter  Estimation 

A  number  of  problems  are  associated  with  recursively  updating  all 
secondary  parameters  according  to  eq.  (3-11)  at  every  sampling  period. 
As  the  number  of  inputs,  outputs,  and  the  order  of  the  linear 
difference  equation  model  increased,  the  computational  effort 
increases.  Any  more  than  two  inputs  and  two  outputs  is  impractical  due 
to  long  convergence  time  and  computation  time.  Also,  the  more 
parameters  which  must  be  updated,  the  more  excitation  the  system 
requires  for  parameter  convergence.  Another  problem  concerns 
convergence  to  the  save  operating  point.  Even  though  each  estimator 
may  begin  at  a  different  operating  point  set  from  a-priori  information, 
the  estimators  will  eventually  converge  to  the  same  model  [13]  which 
makes  the  MMA  redundant. 

IV. A. 2.  Partial  Secondary  Parameter  Estimation 

The  secondary  estimators  may  be  partially  estimated  as  explained  in 
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Section  2.  This  method  requires  much  less  computational  effort  than 
full-scale  parameter  estimation.  Also,  since  fewer  parameters  must  be 
estimated,  they  should  converge  faster.  The  common  convergence  problem 
of  each  estimator  mentioned  with  full-scale  estimation  is  also  avoided 
since  some  parameters  remained  fixed. 

TV. A. 3.  Fixed  Secondary  Parameter  Estimation 

The  designer  also  has  the  option  of  leaving  all  parameters  fixed  at 
a  nominal  flight  condition,  skipping  secondary  parameter  estimation 
altogether.  This  method  is  very  simple,  requiring  relatively  little 
computational  effort.  Its  only  major  drawback  is  that  it  is  only  able 
to  estimate  parameters  at  discrete  flight  conditions.  This  is 
acceptable,  however,  given  a  robust  control  law. 

IV. B.  MMA  Implementation 

Output  coupling  effects  upon  parameter  estimation  were  neglected  in 
t  all  simulations.  Therefore  the  MMA  weighting  coefficient  update 

,  eq.  (3-12)  is  modified  as  a  scalar  update  equation: 

f 

t 


ai(k)  =  CsajJk  -  1)  [diag{Vi(k)  )]-V2exp{-l/2SUM[e2ifj(k) 


IV. C.  Weighting  Coefficient  and  Prediction  Error  Variance  Estimation 
Inspection  of  eq.  (3-12)  shews  that  the  weighting  coefficient 
update  is  a  Gaussian  distribution  where: 

aiM/Cgaifk  -  1)  =  p(ei)  =  *N(0,Vi(k))  (4-4) 

in  where 

vi(k)  =  v(k)  +  oT(k)Pi(k)o(k)  (4-5) 

Estimates  of  v^(k)  for  use  in  eq.  (3-10)  through  eq.  (3-12)  can  be 
based  on  either  e^(k)  or  Pf(k),  for  example 


k 

Viest(k)  =  (i/k)SUM{ei(ii)eiT(ii))  (4-6) 

ii=l 

Viest(k)  =  v(k)  +  e>T(k)pi(k)c3(k)  (4-7) 

where 

k 

v(k)  =  (l/k)SUM{ej (ii)ejT(ii) )  -  oT(k) Pj (k) o(k)  (4-8) 

ii=l 

and  j  is  the  argument  of  the  largest  a^(k)  and  v^est(k)  is  an  estimate 


of  v^(k).  Eq.  (4-6)  is  based  on  e^(k)  and  is  similar  to  the  estimate 
of  v(k)  in  Section  2.  Eq.  (4-7)  requires  a  noise  variance  term,  v(k) , 
which  is  calculated  in  eq.  (4-8)  and  is  the  same  for  each  weighting 


.  *  .  - 
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coefficient.  Differences  in  Vi(k)  are  based  on  differences  in  Pi(k). 

A  physical  interpretation  of  the  two  terms  in  eq.  (4-7)  is  that  of 
a  separation  of  modeling  and  noise  errors.  Define  the  parameter  error 
as 


Wi(k)  =  ©i(k)  -  0(k)  (4-9) 

then 

ei(k)  =  y(k)  -  ©T(k)e(k)  -  oT(k)w^(k)  (4-10) 

The  first  two  terms  represent  noise  due  to  actuators,  sensors  and 
linear  approximation  errors.  The  third  term  is  an  additional  model 
error  term  induced  by  the  model  not  being  the  best  fitting  model. 
Therefore  when  neglecting  cross-correlation  terms  and  dropping  the  time 
index,  k 


E{ejeiT)  =  E([y  -  0T©] [y  -  0T©]T)  +  E{ [oTwi] [oTwi]T)  (4-11) 

=  E{[y  -  0T©][y  -  0T© ] T }  +  E{0TwiWiT0)  (4-12) 

=  E{  [y  -  0T©]  [y  -  0T©]T)  +  0TE{w-jw;jT}0  (4-13) 

=  V(k)  +  0TPj0  (4-14) 


An  alternative  to  eqs.  (4-6)  and  (4-7)  is  to  implement  a  constant 
v^(k)  in  eqs.  (3-10)  -  (3-12).  Andersson  [2]  claims  the  MMA  is 
relatively  insensitive  to  vj_es'*:  (k)  and  a  close  guess  may  be  good  enough 
when  v(k)  is  constant.  Andersson1 s  simulations  show  acceptable  results 
for  an  assumed  value  of  v(k) 

Though  commonly  used  in  literature,  "noise  variance"  is  a 
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misleading  name  for  v(k)  since  v(k)  is  composed  not  only  of  sensor 
noise  variance  but  also  of  linear  approximation  errors.  In  practical 
situations,  sensor  noise  will  be  much  larger  than  linear  approximation 
errors  and  the  linear  approximation  errors  are  neglected.  In  computer 
simulations  which  assume  no  sensor  noise,  however,  linear 
approximations  errors  can  not  be  neglected. 

IV. D.  Improving  Vj^^fk)  Via  Input/Output  Filtering 

Input-output  filtering  is  necessary  in  adaptive  control 
applications  in  order  to  estimate  and  track  v^(k) .  Filters  which 
estimate  v^(k)  may  be  limited  by  the  assumption  that  vj(k)  is 
relatively  slowly  time  varying.  Low-pass  filtering  has  the  effect  of 
reducing  the  variation  in  time  of  v(k) .  Since  v^(k)  is  a  function  of 
0(k),  reducing  the  high  frequency  variation  in  u(k)  and  y(k)  will 
reduce  the  variation  of  v-jjk)  .  A  high-pass  filter  may  be  required  to 
remove  parameter  estimate  bias  when  u(k)  is  non-zero  mean  [10,20]. 
Although  high-pass  filtering  may  not  be  useful  in  conditioning  v^(k) , 
its  possible  bias  on  vjjk)  should  be  considered  upon  its 
implementat ion . 

IV. E.  Convergence  Analysis 

Although  it  is  difficult  to  make  any  general  conclusions  about 
parameter  estimate  convergence  time  given  certain  assumptions  a 
convergence  analysis  may  be  possible  and  useful  when  analysing 
simulation  data.  For  example,  a  convergence  analysis  of  the  MMA  tor 
the  assumption  that  the  prediction  error  vari ance  estimates  are  to 


large  and  time  invariant  may  provide  information  relating  to  the 
sensitivity  of  the  algorithm  to  the  value  of  the  estimated  prediction 


error  variance. 

Example  4-1 

This  example  calculates  an  approximate  convergence  time  for  a  SISO 
adaptive  simulation  which  uses  a  2-model  MMA  estimator.  Assume  the  MMA 
weighting  coefficients  have  reached  steady  state  values  a^(k-n)  =0.01 
and  a2(k-n)  =0.99.  A  maneuver  is  then  simulated  by  changing  the  plant 
dynamics  such  that  at  time  kT  the  weighting  coefficients  have  changed 
to  a^(k)  =  0.99  and  a2(k)  =0.01.  In  such  case,  the  MMA  weighting 
coefficient  update  equations  are 

ai(k)  =  [ai(k  -  l)/[a1(k)  +  a2(k)]]  v^V^expie]2 (k)/v2 (k)  )  (4-15a) 
a2(k)  =  [a2(k  -  l)/[ai(k)  +  a2(k)]]  v2~1/2exp{e12 (k)/v2 (k)  }  (4-15b) 

If  E{e^2(k) >  »  v^(k)  and  E{e22(k) }  »  v2 (k)  then  eq.  (4-15)  simplies 
such  that 


ai(k)  =  [a1(k  -  l)/[a!(k)  +  a2(k)]]  v1“V2 
a2 (k)=[a2 (k  -  l)/[a1(k)  +  a2(k)]]  v2_1/2 


and 


a! (k)  ax(k  -  l) v2 

al (k  -  n) 

v2n/2 

a2  (k)  a2  (k  -  ljvj1/2 

a2 (k  -  n) 

v  ]  h/2 

(4-16a) 

(4-16b) 


(4-17) 


n/2 


log[a1(k)a2(k  -  nj/^kja^k  -  n)  ] 


(4-18) 


log[v2(k)/v1(k) 

IV. F.  System  Filters 

Three  filters  are  incorporated  into  the  system  shewn  in  Figure  (IV- 
1)  to  smooth  variations  in  the  parameter  estimates.  The  first  is  a 
digital  band-pass  filter  which  filters  y(k)  and  u(k) .  The  non-linear 
second  filter  is  added  to  the  weighting  coefficients  within  the  MMA. 

It  smooths  sudden  changes  in  the  model  probabilities  which  cause 
destablizing  rapid  changes  in  the  final  parameter  estimates.  The  third 
filter  smooths  estimates  before  entering  the  control  law  algorithm. 

IV.F.l.  Input-Output  Filter 

The  input-output  signals  are  filtered  by  a  band-pass  filter.  The 
low  frequency  components  of  these  signals  must  be  removed  to  reduce 
parameter  estimate  bias,  while  high  frequency  signal  components  are 
removed  to  smooth  input  excitation  [6,10,11,20].  The  filter  is 
mechanized  as  a  sixth-order  butterworth  digital  band-pass  filter. 

IV. F. 2.  Weighting  Coefficient  Filter 

The  weighting  coefficients,  ai(k),  are  filtered  to  limit  the  rate 
at  which  they  can  change  in  a  given  sampling  period.  This  corresponds 
to  a  limit  of  the  amount  of  old  information  which  is  thrown  away  from 
ai(k)  during  an  update.  Large  variations  in  a [ (k)  in  turn  cause  large 
variations  in  the  parameter  estimates  which  is  destabilizing. 


IV. F. 3.  Rate  Limiting  Filter  Design 


A  rate  limiting  filter,  which  satisfies  the  listed  design  r  :t.  :  .  i 
and  requires  no  renormalization,  limits  the  change  from  a i ( k  -  l,  t.. 
a£(k)  such  that: 

ai,fil(k)  =  ai(*  -  1)  +  c[dai(k)]  (4-19) 

where 


dai  =  ai(k) 


a-p(k  -  1)  and 


c  = 


1 

amax/aj 


if  %ax  ^  aj 
if  amax  -  aj 


and 

j  =  arguement  of  the  maximum  da^ 

a^ax  =  a  constant  representing  the  largest  value  of  dag (k)  permitted. 


Summary 

Implementing  the  MMA  requires  develcping  a  specific  mechanization 
and  choosing  appropriate  filters.  The  MMA  is  very  flexible  in  that  it 
contains  many  application  dependent  methods  of  calculating  key  vari¬ 
ables  such  as  secondary  parameter  estimates  and  prediction  error  vari¬ 
ances.  The  following  section  will  experimentally  examine  many  of  these 


issues . 


The  objective  of  this  section  is  to  experimentally  determine  which 
of  the  mechanizations  developed  in  Section  4  works  best  for  in-flight 
simulation  by  judging  them  on  their  effect  on  three  main  design  crite¬ 


ria  for  the  MMA.  Convergence  to  the  best  fitting  model  is  the  first 
objective.  Even  though  the  control  law  is  robust,  the  controller's 
tracking  ability  will  degrade  given  biased  estimates.  Fast  parameter 
tracking  is  the  second  objective.  Slew  parameter  changes  which  occur 
when  the  aircraft's  flight  condition  changes  must  be  followed  by  the 
estimator.  The  third  and  final  objective  is  for  the  parameter 
estimates  to  be  insensitive  to  sensor  noise.  In  other  words,  the  first 
two  objectives  hold  when  a  realistic  amount  of  sensor  noise  is  added. 

V.B.  Computer  Implementation 

Simulations  are  run  on  a  VAX  11/780  digital  computer  using  Fortran 
subroutines  interfaced  with  the  Matrixx  computer-aided  control  design 
software  package  [ 9 ] . 

V.C.  Set-Up 

The  simulations  approximate  the  motion  of  an  *AFTI  F-16  aircraft 
by  linearized  longitudinal  equations  of  motion  given  in  state  space 


*  The  VISTA,  currently  under  development,  will  be  represented  by  the  AFTI  F- 
16  which  has  many  of  the  control  performance  features  of  the  VISTA.  The  AFTI 
F-16  is  discussed  in  appendix. 
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form  in  appendix  B.  Assume  at  time  t=0  the  host  is  tracking  * itself 
at  0.60  Mach,  10,000  ft.  MSL.  The  reference  tracking  signals  are  from 
a  real-time,  non-linear  simulation  [19,20]  with  elevators  and  flaperons 
as  control  surfaces  with  flight  path  angle  and  pitch  rate  as  outputs. 

At  time  t  =  6  seconds  the  linearized  equations  are  changed  to 
approximate  0.31  Mach,  10,000  ft.  MSL  to  represent  a  worst  case 
parameter  change.  The  MMA  uses  2  estimators.  The  linear  difference 
equation  parameters  of  estimator  1  are  initialized  at  0.60  Mach,  10,000 
ft.  MSL  (appendix  B) .  The  parameters  of  estimator  2  are  initialized  at 
0.31  Mach,  10,000  ft.  MSL.  The  weighting  coefficients,  ajjk), 
initially  are  each  0.5.  These  quantities  are  fixed  until  the  algorithm 
is  turned  on  at  t  =  2  seconds. 

V.D.  Sensor  Noise 

Sensor  noise  is  considered  in  some  simulations.  The  form  and 
magnitude  of  the  noise  is  discussed  in  appendix  B. 

V.E.  Tuning  Haqqlund ' s  Algorithm 

Since  fault  detection  is  not  being  used  the  asymptotic  covariance 
variable,  a,  is  the  only  important  tuning  parameter  which  needs 
discussion.  Hagglund's  algorithm  (appendix  C)  seeks  to  force  the 
parameter  covariance  matrix,  Pi(k),  to  converge  to  "al",  where  "I"  is 
the  identity  matrix,  when  parameters  remain  time  invariant.  Therefore 
since  an  initial  estimate  of  the  parameter  vector  is  available  a- 

*  The  model  aircraft  in  each  simulation  was  the  same  as  the  host  aircraft  to 
avoid  implementing  model-follcwing  equations.  Model-following  theory  for  in¬ 
flight  simulation  is  well  established  in  the  literature  [20]. 
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priori,  a  value  of  "a"  should  be  based  upon  it.  Note  also  since  "a" 
directly  effects  the  asymptotic  parameter  covariance  it  can  be  used  to 
adjust  the  sensitivity  of  the  estimator.  Reducing  the  value  of  "a" 
reduces  the  parameter  estimate  noise  while  increasing  "a"  allows  the 
estimator  to  respond  more  quickly  to  changes  in  aircraft  dynamics. 

V.F.  Parameter  Estimator  Variations 

Table  (V-l)  lists  ten  possible  variations  of  the  estimator  shown  in 
Figure  (IV-1) .  For  the  advantages  shewn,  this  experiment  is  limited  to 
method  1  of  estimating  the  prediction  error  variance  in  conjunction 
with  partial  secondary  estimation  with  Hagglund's  algorithm  and  no 
secondary  parameter  estimation. 


Table  (V-l)  comparison  of  Secondary  Parameter  Estimation  Methods 


Secondary  Method  1  for 
Estimation  Estimating  Vj(k) 

Full  Scale  -High  Computation  Cost 
RLS  -Parallel  Estimators  Converge 

to  Same  Operating  Point 


Method  2  for 
Estimating  v^fk) 

-High  Computation  Cost 
-Parallel  Estimators  Converge 
to  Same  Operating  Point 


Full  Scale  -High  Computation  Cost 
Hagglund  -Impractical  for  High  Order 

-Parallel  Estimators  Converge 
to  Same  Operating  Point 


-High  Computation  Cost 
-Same  as  Method  1  and 
-Modifications  to  P  destroy 
goodness  of  fit  information 


Partial 

RLS 


-Acceptable  Computation  Cost  -P  only  Partially  Estimated 
Continuous  Envelope  Coverage 


Partial  -Same  as  Partial  RLS  but 
Hagglund  better  Wind-up  Resistance 


-Same  as  Partial  RLS  but 
better  Wind-up  Resistance 


No  Update  -Small  Computation  Costs 

-Discrete  Operating  points 


-Not  Possible, 

P  not  estimated 
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V.G.  Estimator  Set-up 

The  experiment  consists  of  four  simulations  of  the  system,  each 
with  a  different  method  of  parameter  estimation: 

Simulation  1  -  Partial  secondary  parameter  estimation  via  Hagglund's 
algorithm,  error  variance  estimation  via  method  1,  no  noise. 

Simulation  2  -  Partial  secondary  parameter  estimation  via  Hagglund's 
algorithm,  error  variance  estimation  via  method  1,  with  noise. 

Simulation  3  -  No  secondary  parameter  estimation,  error  variance  esti 
mation  via  method  1,  no  noise 

Simulation  4  -  No  secondary  parameter  estimation,  error  variance  esti 
mation  via  method  1,  with  noise. 

V.H.  Multiple  Model  Estimator  Performance 

Simulation  data  shows  the  MMA  is  useful  for  adaptive  control  for 
in-flight  simulation  by  meeting  the  stated  objectives.  Plots  of  param¬ 
eter  estimates  show  parameter  convergence  and  tracking  after  a  jump 
change  in  host  aircraft  dynamics.  Output  data  plots  show  almost  per¬ 
fect  tracking  of  the  model  aircraft  output  by  the  host  aircraft. 


V. I .  Comparative  Estimator  Variation  Performance 

Examination  of  experimental  data  provides  several  important  guide¬ 
lines  for  using  MMA  variations  shewn  in  Table  (V-l) . 

1)  Prediction  error  variance  estimation  method  1  (section  3)  is 
preferable  to  method  2  in  most  instances.  This  is  especially  true 
if  Hagglund's  algorithm  is  used  for  secondary  parameter  estimation. 

2)  Secondary  parameter  estimates  don't  converge  in  general  to  the 
best  fitting  model  when  partial  parameter  estimation  is  used. 

3)  The  key  to  MMA  performance  is  the  prediction  error  variance 
estimates,  v^(k).  If  prediction  error  distributions  of  each  model 
have  similar  variances,  the  weighting  coefficients,  a^(k),  may 
converge  to  the  wrong  model. 
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4)  Hie  MMA's  tolerance  to  sensor  noise  is  dependent  on  input  excita¬ 
tion.  Sensor  noise  variance  must  remain  smaller  than  model  error 
variance  to  ensure  probability  convergence. 

5)  Proper  input/output  filtering  is  essential. 

6)  Weighting  coefficients,  ap(k),  tend  to  converge  to  unity  for 
which  ever  model  has  the  smallest  prediction  error  variance  esti¬ 
mates. 

7)  Input  excitation  has  a  direct  effect  on  the  prediction  error 
variance  due  to  model  error. 

Guideline  1 

Plots  of  elements  of  Pp(k)  illustrate  the  problem  associated  with 
simultaneously  using  method  2  to  estimate  vp(k)  when  using  Hagglund's 
algorithm  for  secondary  parameter  estimation.  Notice  the  plots  of 
element  (1,1)  of  various  Pp(k)  matrices,  which  illustrate  that  in 
general,  the  Pp(k)  matrix  whose  arguement  corresponds  to  the  best 
fitting  model  has  the  smallest  elements.  Hagglund's  algorithm,  howev¬ 
er,  forces  the  diagonal  elements  of  Pp(k)  to  converge  to  a  predeter¬ 
mined  constant  (see  discussion  of  Hagglund's  Algorithm  parameter  "a"  in 
Appendix  C) .  If  the  diagonal  elements  of  the  estimated  Pp(k)  matrix 
converge  to  this  constant  in  both  models,  the  portion  of  the  prediction 
error  variance  due  to  modeling  errors  will  be  biased,  degrading  the 
performance  of  the  MMA. 

Guideline  2 

When  using  partial  parameter  estimation,  secondary  estimates  of  Bp 
do  not  converge  to  the  Bp  associated  with  the  best  fitting  model. 
Instead,  they  tend  to  vary  around  the  Bp  associated  with  same  model  as 
the  fixed  parameters.  This  behavior  is  illustrated  in  plots  of 


secondary  B]_  estimates  from  Simulations  1  and  2. 

The  bias  in  means  parameter  estimates  will  be  limited  to  points 
on  the  flight  envelope  for  which  the  secondary  models  are  calculated. 

A  robust  control  law,  however,  should  still  perform  adequately  provided 
the  models  are  not  spread  too  far  apart  within  the  flight  envelope. 
Model  placement  is  an  important  issue  and  is  recommended  as  an  area  of 
continued  research. 

Guideline  3 

The  weighting  coefficients  of  the  MMA  converge  to  the  model  with 
the  smallest  prediction  error  variance  estimate.  Convergence  behavior 
is  difficult  to  determine  when  prediction  variance  estimates  are  of 
similar  magnitude.  Notice  the  plots  of  the  prediction  error  variance 
in  Simulation  2  and  4.  Around  7  seconds,  when  little  input  excitation 
occurs  and  the  variance  of  each  model  is  similar,  the  probability 
calculation  does  not  converge. 

Guideline  4 

The  prediction  error  variance  due  to  modelling  error  must  be  larger 
than  that  due  to  sensor  noise  in  order  for  the  weighting  coefficients 
calculations  to  converge  to  the  best  fitting  model.  If  the  prediction 
errors  are  due  more  to  sensor  noise  than  modeling  error,  then  the 
estimated  prediction  error  variance  will  be  almost  the  same  for  each 
model  .  Then,  from  guideline  3,  probability  calculations  are  suspect. 
This  behavior  is  illustrated  in  prediction  error  variance  estimates  in 


Simulation  2  and  4. 
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Guideline  5 

Proper  input/output  filtering  is  essential  to  proper  MMA  operation. 
Prediction  error  noise  spikes  may  cause  spikes  in  prediction  error 
variance  estimates  which  may  cause  weighting  coefficients  to  diverge. 
This  behavior  is  observed  in  simulations  (not  shewn)  where  input  and 
output  measurements  are  not  filtered.  Heavy  filtering,  however,  dis¬ 
cards  good  information  by  reducing  the  difference  between  filtered 
prediction  error  variances.  Filtering  also  causes  a  "coloring"  of  Jie 
prediction  errors  [15].  No  quantitative  method  for  choosing  an  input- 
output  filter  has  been  developed  so  trial  and  error  is  used  to  choose 
the  band  width  of  the  bard-pass  filter  used  in  these  simulations. 

Guideline  6 

Weighting  coefficients  tend  to  converge  very  quickly,  with  the 
weighting  coefficient  of  the  best  fitting  model  converging  to  unity. 
This  behavior  is  demonstrated  in  simulations  (not  shown)  where  the 
actual  operating  point  is  between  the  two  points  used  to  calculate 
initial  secondary  parameter  estimates. 

Guideline  7 

The  theoretical  dependence  of  the  input  on  the  prediction  error 
variance  due  to  modeling  error  is  discussed  in  Section  4.  Comparison 
of  the  input  signal  and  the  prediction  error  variance  plots  verifies 
this  relationship. 
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Table  (V— 2)  I list  of  Plots  -  Experimental  Results 


UPPER 

Secondary  estimated 
B]_,  element  1,1  -  Sim.  1 

Primary  estimated 

Bj,  element  1,1  -  Sim.  1 

Primary  estimated 

B]_,  element  1,1  -  Sim.  3 

Primary  estimated 

Bi,  element  2,2  -  Sim.  1 

Primary  estimated 


Bl, 

Pi, 

element 

element 

2,2 

1,1 

-  Sim.  3 

-  Sim.  1 

Pi, 

element 

2,2 

-  Sim.  1 

vi, 

element 

1,1 

-  Sim.  1 

vi, 

element 

1,1 

-  Sim.  3 

vi, 

element 

2,2 

-  Sim.  1 

Vi, 

element 

2,2 

-  Sim.  3 

al, 

Sim.  1 

al, 

Sim.  3 

Flight  path  angle 
Commanded  and  actual 
Sim.  1 


Flight  path  angle 
Commanded  and  actual 
Sim.  3 

Pitch  rate 
Commanded  and  actual 
Sim.  1 

Pitch  Rate 
Commanded  and  actual 
Sim.  3 


IPWER 

Secondary  estimated 
B]_,  element  1,1  -  Sim.  2 

Primary  estimated 

B]_,  element  1,1  -  Sim.  2 

Primary  estimated 

B^,  element  1,1  -  Sim.  4 

Primary  estimated 

B^,  element  2,2  -  Sim.  2 

Primary  estimated 

B]_,  element  2,2  -  Sim.  4 

Pi,  element  1,1  -  Sim.  2 

Pj,  element  2,2  -  Sim.  2 

V|,  element  1,1  -  Sim.  2 

v-[,  element  1,1  -  Sim.  4 

v^,  element  2,2  -  Sum.  2 

Vj^ ,  el  ement  2,2  -  Sim.  4 

a^,  Sim.  2 

a^,  Sim.  4 

Flight  path  angle 
Commanded  and  actual 
Sim.  2 

Flight  path  angle 
Commanded  and  actual 
Sim.  4 

Pitch  rate 
Commanded  and  actual 
Sim.  2 

Pitch  Rate 

Commanded  and  actual 
Sim.  4 


VI.  Conclusion 


The  purpose  of  this  thesis  is  to  use  the  MMA  for  parameter  identi¬ 
fication  used  in  adaptive  control  of  an  in-flight  simulator.  The 
thesis  is  an  extension  of  Pineiro's  thesis  [20]. 

Pineiro  sought  to  control  an  in-flight  simulator  via  an  adaptive 
model-following  PI  control  law  [20-24].  The  control  law  bases  its 
control  gains  upon  the  parameters  of  a  linear  difference  equation  model 
which  describes  the  input/output  relationship  of  the  host  aircraft's 
dynamics.  The  parameters  of  this  model  are  estimated  from  input  and 
output  measurements  via  a  recursive  parameter  estimation  algorithm. 

This  thesis  differed  from  Pineiro's  in  that  the  recursive  parameter 
estimation  algorithm  is  replaced  by  the  MMA  to  increase  convergence 
speed  to  the  best  fitting  model  parameters.  Slow  convergence  is  a 
critical  problem  for  recursive  parameter  estimation  algorithms  when 
used  with  higher  order  systems  such  as  those  typically  required  to 
describe  aircraft  dynamics. 

A  computer  simulation  is  performed  which  shows  excellent  results. 
The  MMA  converge  quickly  to  the  best  fitting  model.  The  host  aircraft 
tracks  the  model  aircraft  accurately.  Degraded  but  acceptable  perform¬ 
ance  results  when  sensor  noise  was  added. 

Experimental  data  also  reveals  useful  information  on  tuning  the  MMA 
for  actual  use.  The  MMA  is  flexible  and  can  be  implemented  in  several 
ways.  Implementation  options  are  discussed  and  performance  compared. 
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Recommendations  for  Further  Study 

Optimal  placement  within  the  flight  envelope  is  one  example  re¬ 
search  area  which  needs  to  be  investigated.  The  computational  cost  of 
adding  more  a-priori  models  to  the  simulation  is  a  trade-off  of  more 
accurate  parameter  estimation.  A  moving  bank  scheme  [16]  which  varies 
the  models  used  in  the  MMA  may  be  used  to  reduce  the  number  of  models 
needed  to  cover  the  entire  envelope. 

More  research  also  is  required  to  decrease  the  MMA's  sensitivity  to 
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sensor  noise.  One  possibility  is  to  turn  the  MMA  off  during  noisy 
periods  and  when  the  system  is  receiving  little  input  excitation. 


Figure  (A-l)  AFTI  F-16 

The  AFTI  (Advanced  Fighter  Technology  Integration) /F-16  aircraft, 
shown  in  Figure  (A-l) ,  is  an  F-16A  air  superiority  fighter  modified  to 
be  a  test  bed  for  evaluating  new  aircraft  technologies.  These 
modifications  include  the  addition  of  two  vertical  canards  which  are 
mounted  on  the  engine  inlet,  control  surfaces  that  allow  independent 
motion  of  the  trailing  edge  flaps  and  independent  motion  of  the 
horizontal  tail  halves,  and  a  redundant,  digital  fly-by-wire  flight 


control  system. 


Primary  mission  tasks  involve  flight  maneuvers  which  are  divided 
into  four  flight  control  configuration  modes.  The  first  is  referred  to 
as  the  normal  mode  and  is  used  for  takeoff,  cruise,  landing,  and  air 
refeuling  tasks.  The  second  is  the  air-to-air  gunnery  mode,  used  for 
precise  target  pointing  and  evasive  maneuvering.  The  third  is  the  air- 
to-surface  gunnery  mode,  used  for  precise  pointing  needed  for  strafing 
ground  targets.  The  final  configuration  is  the  air-to-surface  bombing 
mode  used  for  precise  velocity  control  required  for  accurate  bombing. 
The  aircraft  models  used  in  this  thesis  were  obtained  when  the  aircraft 
was  in  the  normal  mode. 


1  Data 


These  are  the  continuous  time  state  space  model  approximating  the 
equations  of  motion  used  in  the  simulation: 


=  [U!(t)  u2(t)]T 
=  [Yi(t)  Y2(t)]T 


u3(t.)  =  elevator  defection  (degrees) 
u2(t)  =  tlaperon  deflection  (degrees) 


y3(t)  =  flight  path  angle  (degrees) 
y2(t)  =  pitch  rate  (degrees/sec) 


=  [Xx(t)  x2(t)  x3(t)  x4(t)]T 


x3(t)  =  pitch  angle  (degrees) 
x2 (t)  =  perturbation  velocity 
(ft/sec) 

x3(t)  =  angle  of  attack  (degrees) 
x4 (t)  =  pitch  rate  (deqrees/sec) 
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0  0  0  1 


=  output  matrix 


Nominal  Flight  Condition  =  10,000  ft.,  MACH  .31 


0.0000 
-3 . 1884D+1 
A  =  -1.2770D-2 

7 . 5310D-4 


0.00000 
1.66435 
-7 . 70780D-2 
-3.25199 


0.0000 
-1.2236D-2 
-2 . 9100D-4 
6. 0000D-5 


0.00000 
-4.29844 
-6.91360D-2 
3 . 25307D-1 


0.0000 
1 . 7789LH-1 
—4 . 8935D-1 
1.78948 


1.0000 
-4 . 5361D+1 
0.999914 
-3 . 8710D-1 


Nominal  Flight  Condition  =  10,000  ft.,  MACH  .60 


0.0000 
-3 . 2046D+1 
A  =  -7 . 0900D-3 

5.8700D-4 


0.0000 
-1.5269D-2 
-2 . 0800D-4 
1 . 3220D-4 


0.0000 
2 . 3608D+1 
-6 . 4800D-1 
2.28298 


1.0000 
-3.6205CH-1 
0.999924 
-4 . 1422D-1 


.  \  N  *.  \  \  *.  %  *.  %  *,  *. 

.  .  -**»  ■  .  «  ..  -  _  -* 
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0.00000 
2.01664 
-9 . 31110D-2 
-4.72071 


0 . 00000 
-1.30775 
-1.37395D-1 
-2 . 303800-1 


Model  Structure 


y (k)  =  -Axy (k  -  1)  -  A2y(k-2)  -  A3y(k-3)  -  A4y(k-4) 

+  B1u(k  -  1)  +  B2U(k-2)  +  B3u(k-3)  +  B4u(k-4)  +  e(k) 


Parameters  at  Nominal  Flight  Condition:  .60  MACH;  10,000  ft.  MSL 


Element 

Element 

Element 

Element 

Matrix 

llxll 

(1,21 

(?,11 

(.2,2) 

B1 

1.411739D-3 

2 . 386221D-3 

-1.324101D-1 

-1.311730D-2 

A1 

-3.982372 

0.0 

0.0 

-3.982372 

B2 

-4 . 241853D-3 

-7 . 145616D-3 

3 . 956153D-1 

3 . 907266D-2 

A2 

5.946657 

0.0 

0.0 

5.946657 

B3 

4.232493D-3 

7.129810D-3 

-3 . 940006D-1 

-3 . 879344D-2 

A3 

-3.946199 

0.0 

0.0 

-3.946199 

B4 

-1.402381D-3 

-2.370406D-3 

1 . 307954D-1 

1.283808D-2 

A4 

9 . 819136D-1 

0.0 

0.0 

9.819136D-1 

Parameters  at  Nominal  Flight  Condition:  .31  MACH;  10,000  ft.  MSL 


Matrix 

Element 

(l.D 

Element 

1.1,21 

Element 

12x11 

Element 

(2,2) 

B1 

1. 000519D-3 

1. 504922D-3 

-5. 588077D-2 

-2 . 567509D-3 

A1 

-3.988486 

0.0 

0.0 

-3.988486 

B2 

-3 . 001350D-3 

-4 . 508592D-3 

1. 672050D-1 

7.644987D-3 

A2 

5.965235 

0.0 

0.0 

5.965235 

B3 

2.996832D-3 

4 . 501846D-3 

-1 . 667678D-1 

-7 . 587455D-3 

A3 

-3.965013 

0.0 

0.0 

-3.965013 

B4 

-9 . 960008D-4 

-1.498176D-3 

5 . 544359D-2 

2 . 509977D-3 

A4 

9 . 882635D-1 

0.0 

0.0 

9 . 882635D-1 

Noise 

Sensor  noise  was  generated  by  a  random  number  generator,  where  the 
random  number,  r^,  fit  the  bound  0  <  rj_  <  1.  In  other  words,  p(r^) 
*U(0, 1) ,  where  E{rjJ  =  .5,  E{ri2}  =  1/3,  Var{rjJ  =  1/12  [17]. 


A  pdf  with  statistics  similar  to  a  gaussian  distribution  can  be 
generated  as 

12 

W  =  SUM(ri)  -  6 

j=l 

where  E{W)  =  0,  E{W2}  =  Var{W}  =  1 

The  variance  of  W  was  calculated  from  the  following  theorem  [17]: 

If  {r]_,  r2,...,ri}  is  a  sequence  of  independent  random  variables,  then 

i  i 

Var(SUM{rj  })  =  SUM{Var(rj)  } 

j=l  j=l 

Simulated  sensor  noise  is  made  by  scaling  W  with  the  desired  noise 
variance. 

Noise  Figures 

Zero-mean  pseudo  gaussian  noise  of  the  following  strength  was  added 
to  the  state  vector: 


Variance  State 

7.438X10-11  x1 

0.000  X2 

4.896X10-10  x3 

3.404X10-9  X4 


*  p(x)  =  U(L,U)  is  a  uniform  distribution,  where  p(x)  =  0  for  x  •  I., 
and  p(x)  =  [U-L]-1  for  L  <  x  <  U. 


ndix  C  -  Haqqlund's  Alqorithrti  for  Tracki 


Changes 


The  goal  of  the  estimator  developed  by  Hagglund  is  to  weight  the 
incoming  information  so  that  the  covariance  matrix  becomes  preport ional 
to  the  identity  matrix.  The  diagonal  elements  of  P(k)  may  be 
interpreted  as  approximations  of  the  variances  of  the  corresponding 
parameters.  Therefore,  P-1(k)  is  a  measure  of  "goodness  of  fit" 
which  should  be  kept  constant.  In  other  words,  if  no  information  is 
in  coming  nothing  is  forgotten  but  if  the  information  content  of  the 
current  measurement  is  large,  old  information  is  discarded  quickly  for 
fast  parameter  adaptation. 

If  information  is  to  be  discounted  according  to  the  new  principle 
the  P(k)  update  equation  must  be  modified.  The  RLS  P(k)  update,  eq. 
(2-31) ,  can  be  alternately  expressed  as  a  P-1(k)  update 


P-1  (k)  =  P^-(k  -  1)  +  v_1(k)0(k)0T(k) 


(C-l) 


Eq.  (C-l)  discounts  old  information  exponentially  in  all  directions. 
Hagglund' s  algorithm  modifies  eq.  (C-l)  as 

P_1(k)  =  P-1(k  -  1)  +  v_1(k)0(k)0T(k)  -  a(k)0(k)oT(k)  (C-2) 

=  P-1(k  -  1)  +  [v-1(k)  -  a(k) ]0(k)oT(k)  (C-3) 

■where  a(k)  is  a  discounting  factor.  The  new  information  is 
proportional  to  0(k)o^(k)  and  it  may  be  said  that  the  new  informatior 


is  coming  in  the  direction  of  0(k)  while  old  information  is  discounted 
in  the  same  direction.  Eq.  (C-2)  can  be  transformed  via  the  matrix 
inversion  lemma  into  a  P(k)  update  equation: 

P(k  -  l)o(k)oT(k)P(k  -  1) 

[V1(k)  -  a(k)  ]-1  +  oT(k)P(k  -  l)o(k) 

Since  the  form  of  the  covariance  matrix  update  has  changed,  the 
equation  for  updating  the  parameter  estimates  will  also  change.  Using 
eq.  (C-4)  along  with  the  basic  definition  of  the  parameter  update 
equation  in  the  least-squares  algorithm,  Hagglund  derives  the  following 
parameter  estimate  update  equation: 

0(k)  =  6(k  -  1)  +  (l/v(k))P(k)o(k)e(k)  (C-5) 

It  remains  to  be  shewn  hew  to  select  an  appropriate  discounting 
factor  a(k).  Eq.  (C-2)  and  eq.  (C-3)  show  that  a(k)  must  be  positive 
or  information  would  be  added  instead  of  removed  but  if  a(k)  is  too 
large  the  covariance  matrix  could  become  non-positive  definite. 

Hagglund  performed  a  stability  investigation  shewing  that  the 
covariance  matrix  remains  positive  definite  if  a(k)  is  chosen  such 
that 


0  <  a(k)  <  [oT(k)P(k  -  l)o(k)]-1 


(C— 6) 


Furthermore,  in  order  to  obtain  a  diagonal  P-Matrix  of  the  form  "al" 
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where  "a"  is  the  desired  variance  of  the  parameter  estimates  and  "I" 
is  the  identity  matrix,  Hagglund  shows  that  a(k)  must  be  selected  so 
that 


oT(k)P(k  -  i)P(k)P(k  -  l)®(k) 

a(k)  - -  (C-7) 

oT(k) P(k  -  1) P(k  -  i)«(k) 


substituting  eq.  (C-4)  into  eq.  (C-7)  yields  the  following  desired 
value  of  a(k) 


delta^(k) 

ad(k)  =  v-1  (k)  + -  (C-8) 

delta^ (k)oT(k)P(k  -  1)0  (k)  -  1 


1 

^  0T(k)p2(k  -  l)0(k) 


Eq.  (C-7)  can  not  always  be  used,  however,  due  the  restrictions 
given  in  eq.  (C-6) .  Hagglund  shews  that  by  incorporating  the  bounds  of 
a(k)  in  eq.  (C-6)  in  conjunction  with  eq.  (C-8)  the  choice  of  a(k) 
becomes 


0 

if 

ad(k) 

< 

0 

(C-10) 

adW 

if 

0 

< 

ad(k) 

<  C 1/n (k) ] 

l/n(k) 

if 

l/n(k) 

< 

^d(i) 

<  v“l(k)  + 

l/n(k) 

0 

if 

ad(k) 

> 

v-1 (k) 

+  1,/n  (k) 

0T(k)P3(k  -  1)0 (k) 
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Abstract 


Adaptive  control  of  aircraft  model-following  systems  has  shewn 
premising  results  for  in-flight  simulation,  but  the  computational 
expense  and  slow  convergence  of  conventional  parameter  estimation 
techniques  for  higher  order  models  inhibits  their  direct  use  for  in¬ 
flight  simulation.  Computer  simulations  of  adaptive  systems  usually 
assume  some  knowledge  of  model  parameters  in  order  to  maintain  tracking 
fidelity  at  a  reasonable  computational  cost  as  parameters  change.  This 
thesis  incorporates  a-priori  information  into  a  multiple-model  estima¬ 
tion  algorithm  which  assigns  a  probability  weighting  of  each  estimator 
within  a  "bank"  of  estimators.  Final  parameter  estimates  used  in 
adaptive  control  are  formed  as  a  probabalistic  weighted  sum  of  individ¬ 
ual  estimates.  Simulations  of  the  system  show  excellent  tracking 
performance  throughout  the  flight  envelope.  A  moving  bank  scheme  for 
use  over  a  wide  range  of  flight  conditions  is  recommended  as  a  further 
area  of  study. 
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